In [31]:
%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
In [32]:
# use our cdist function developed in the previous notebook
def cdist(I, J):
i, j = np.mgrid[:I.shape[0],:J.shape[0]]
return np.sqrt(np.sum((I[i] - J[j])**2,axis=2))
In [33]:
# create some randomly-located points
n_points=10
I, J = np.random.uniform(size=(2,n_points))
IJ = np.vstack((I,J)).T
plt.figure(figsize=(7,7))
plt.gca().set_aspect('equal')
plt.scatter(I,J)
Out[33]:
In [34]:
# compute the pairwise distances between all the points
D = cdist(IJ,IJ)
D
Out[34]:
In [35]:
# if we sort along one axis we get the distance to the nearest neighbors
# of each point
# note that the distances in the first column are zero because each point
# is closest to itself
nn_d = np.sort(D,axis=1)
nn_d
Out[35]:
In [36]:
# to get the closest k distances, we just truncate this array along the second axis
nn_d[:,:3]
Out[36]:
In [37]:
# argsort of the distance array tells us *which* points are nearest.
# note that the first column is sequential because each point is closest to itself
nn_ix = np.argsort(D,axis=1)
nn_ix
Out[37]:
In [38]:
# and we can subset it by whatever k we want
k=3
knn_ix = nn_ix[:,:k]
knn_ix
Out[38]:
In [39]:
# to get the distances that correspond to the indices, we need
# to use a fairly tricky indexing scheme.
In [40]:
# first create a matrix the same size as knn_ix of row numbers
row, _ = np.mgrid[:IJ.shape[0],:k]
row
Out[40]:
In [41]:
# now we can use the row matrix as the row index into D,
# and the knn_ix matrix as the column index
knn_dist = D[row,knn_ix]
knn_dist
Out[41]:
In [42]:
# here is the resulting function, allowing us to specify k
def knn(query_pts, pts, k=3):
# brute force k nearest neighbors
D = cdist(query_pts, pts)
knn_ix = np.argsort(D,axis=1)[:,:k]
row, _ = np.mgrid[:query_pts.shape[0],:k]
knn_dist = D[row, knn_ix]
return knn_dist, knn_ix
In [43]:
from scipy import stats
# now make a grid
resolution = 250
m = np.linspace(0, 1, resolution)
n = np.linspace(0, 1, resolution)
M, N = np.meshgrid(m,n)
MN = np.dstack((M,N)).reshape(-1,2)
# now find the distasnce to the nearest neighbor of each grid point
dist, nabe = knn(MN,IJ,k=1)
# turn that into an image
dist_img = dist.reshape((resolution,resolution))
plt.figure(figsize=(7,7))
plt.gca().set_aspect('equal')
plt.scatter(I, J, c='black')
# we negate dist_img so that low distances are brighter
plt.imshow(-dist_img,origin='lower',cmap='gray',extent=[0, 1, 0, 1])
Out[43]:
In [44]:
# now let's visualize which neighbor is nearest to each point
# reshape the neighbor array into a label image
# a label image has an integer for each pixel identifying its class
# which in this case is the index of the neighboring point in the original
# set of points
label_img = nabe.reshape((resolution,resolution))
label_img
Out[44]:
In [45]:
from skimage.color import label2rgb
# now randomly pick colors for each value in the label image
colors = np.random.uniform(low=0.25,size=(np.max(label_img)+1,3))
plt.figure(figsize=(7,7))
plt.gca().set_aspect('equal')
plt.scatter(I, J, c='black')
# label2rgb turns a label image into an rgb image given colors
plt.imshow(label2rgb(label_img,colors=colors),interpolation='nearest',origin='lower',extent=[0, 1, 0, 1])
Out[45]:
This is a Voronoi diagram